#front matter
rm(list=ls())
library(foreign)
library(MCMCpack)
library(xtable)
library(lattice)
options(scipen=10)

#set working directory
#setwd('/Volumes/MONOGAN/redistributiveSpending/')
#setwd('E:/redistributiveSpending/')
#setwd('C:/Temp')

#open data
fiscal<-read.dta('fiscineq.dta')

###CODE DATA###
fiscal$lnBroad<-log(fiscal$healthhosptotexp_pc+fiscal$publicwelftotalexp_pc)
fiscal$rBroad<-(fiscal$healthhosptotexp_pc+fiscal$publicwelftotalexp_pc)/fiscal$totalexpenditure_pc
fiscal$totdebtoutstand[is.na(fiscal$totdebtoutstand)] <-1
fiscal$lnDebt<-log(fiscal$totdebtoutstand+1)
fiscal$pre10p_p_con<-fiscal$pre10p_p_con/1000
fiscal$pre10p_hh_con<-fiscal$pre10p_hh_con/1000
fiscal$dv.gsp<-fiscal$public_welf_total_exp/fiscal$gspall/1000
fiscal$ratio9050<-fiscal$pre90p_hh_con/fiscal$pre50p_hh_con

###DESCRIPTIVE FIGURE###
#trellis.device("pdf",color=FALSE,file="ratioLine.pdf")
xyplot(rAdjwelftotexptot~year,groups=abb,type='l',data=fiscal,xlab="Year",ylab="Ratio Welfare Spending")
#dev.off()

###DESCRIPTIVE STATISTICS###
#summary(fiscal)
#names(fiscal)
desc.data<-subset(fiscal,select=c(lnWelfTotPC,lnWelfDirPC,rWelfTotExpTot,rWelfDirExpDir,rAdjwelftotexptot,pregini_hh_con,pregini_p_con,gspgrowth,unemploy,lnpop,pop65,nonwhiteprop,gubelection,inst6010_nom,folded_ranney_4yrs,pre10p_hh_con,pre10p_p_con),subset=year>1975&year<2008)
descriptive.output<-cbind(
apply(desc.data,2,mean,na.rm=T),
apply(desc.data,2,sd,na.rm=T),
apply(desc.data,2,min,na.rm=T),
apply(desc.data,2,max,na.rm=T)
)
rownames(descriptive.output)<-c("Log total per capita welfare","Log direct per capita welfare","Ratio total welfare","Ratio direct welfare","Ratio total welfare minus transfers","Household Gini","Personal Gini","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p household income (1000s)","10p personal income (1000s)")
colnames(descriptive.output)<-c("Mean","Std. Dev.","Minimum","Maximum")
xtable(descriptive.output,digits=4,caption="Descriptive Statistics of Variables, 1976-2008",label='descriptives',align="lrrrr")
quantile(fiscal$pregini_hh_con,c(.1,.9),na.rm=T)
#      10%       90% 
#0.4338348 0.5308072 

###BAYESIAN MODELS###
#Current Household Gini#
dep1.hh.mcmc<-MCMCregress(lnWelfTotPC~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep1.hh.mcmc,quantiles=c(.05,.5,.95))
dep2.hh.mcmc<-MCMCregress(lnWelfDirPC~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep2.hh.mcmc,quantiles=c(.05,.5,.95))
#dep3.hh.mcmc<-MCMCregress(lnBroad~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep3.hh.mcmc,quantiles=c(.05,.5,.95))
dep4.hh.mcmc<-MCMCregress(rWelfTotExpTot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep4.hh.mcmc,quantiles=c(.05,.5,.95))
dep5.hh.mcmc<-MCMCregress(rWelfDirExpDir~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep5.hh.mcmc,quantiles=c(.05,.5,.95))
#dep6.hh.mcmc<-MCMCregress(rBroad~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep6.hh.mcmc,quantiles=c(.05,.5,.95))
dep7.hh.mcmc<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.hh.mcmc,quantiles=c(.05,.5,.95))

#Lagged Household Gini#
dep1.hh.lag.mcmc<-MCMCregress(lnWelfTotPC~ lagPregini_hh_con +gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep1.hh.lag.mcmc,quantiles=c(.05,.5,.95))
dep2.hh.lag.mcmc<-MCMCregress(lnWelfDirPC~ lagPregini_hh_con +gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep2.hh.lag.mcmc,quantiles=c(.05,.5,.95))
#dep3.hh.lag.mcmc<-MCMCregress(lnBroad~ lagPregini_hh_con +gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep3.hh.lag.mcmc,quantiles=c(.05,.5,.95))
dep4.hh.lag.mcmc<-MCMCregress(rWelfTotExpTot~lagPregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep4.hh.lag.mcmc,quantiles=c(.05,.5,.95))
dep5.hh.lag.mcmc<-MCMCregress(rWelfDirExpDir~lagPregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep5.hh.lag.mcmc,quantiles=c(.05,.5,.95))#;summary(dep5.hh.lag.mcmc,quantiles=c(.06,.5,.95))$quantiles[2,]
#dep6.hh.lag.mcmc<-MCMCregress(rBroad~ lagPregini_hh_con +gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep6.hh.lag.mcmc,quantiles=c(.05,.5,.95))
dep7.hh.lag.mcmc<-MCMCregress(rAdjwelftotexptot~ lagPregini_hh_con +gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.hh.mcmc,quantiles=c(.05,.5,.95))

#Current Personal Gini#
dep1.p.mcmc<-MCMCregress(lnWelfTotPC~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep1.p.mcmc,quantiles=c(.05,.5,.95))
dep2.p.mcmc<-MCMCregress(lnWelfDirPC~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep2.p.mcmc,quantiles=c(.05,.5,.95))
#dep3.p.mcmc<-MCMCregress(lnBroad~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep3.p.mcmc,quantiles=c(.05,.5,.95))
dep4.p.mcmc<-MCMCregress(rWelfTotExpTot~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep4.p.mcmc,quantiles=c(.05,.5,.95))
dep5.p.mcmc<-MCMCregress(rWelfDirExpDir~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep5.p.mcmc,quantiles=c(.05,.5,.95))
#dep6.p.mcmc<-MCMCregress(rBroad~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep6.p.mcmc,quantiles=c(.05,.5,.95))
dep7.p.mcmc<-MCMCregress(rAdjwelftotexptot~pregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.p.mcmc,quantiles=c(.05,.5,.95))

#Lagged Personal Gini#
dep1.p.lag.mcmc<-MCMCregress(lnWelfTotPC~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep1.p.lag.mcmc,quantiles=c(.05,.5,.95))
dep2.p.lag.mcmc<-MCMCregress(lnWelfDirPC~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep2.p.lag.mcmc,quantiles=c(.05,.5,.95))
#dep3.p.lag.mcmc<-MCMCregress(lnBroad~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep3.p.lag.mcmc,quantiles=c(.05,.5,.95))
dep4.p.lag.mcmc<-MCMCregress(rWelfTotExpTot~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep4.p.lag.mcmc,quantiles=c(.05,.5,.95))
dep5.p.lag.mcmc<-MCMCregress(rWelfDirExpDir~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep5.p.lag.mcmc,quantiles=c(.05,.5,.95))
#dep6.p.lag.mcmc<-MCMCregress(rBroad~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep6.p.lag.mcmc,quantiles=c(.05,.5,.95))
dep7.p.lag.mcmc<-MCMCregress(rAdjwelftotexptot~lagPregini_p_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_p_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.p.lag.mcmc,quantiles=c(.05,.5,.95))
dep1.output.side<-summary(dep1.p.lag.mcmc,quantiles=c(.08,.5,.95))$quantiles[c(2),c(2,1,3)]; dep1.output.side #92%
dep2.output.side<-summary(dep2.p.lag.mcmc,quantiles=c(.14,.5,.95))$quantiles[c(2),c(2,1,3)]; dep2.output.side #86%




#Geweke Diagnostics#
geweke.diag(dep1.hh.mcmc);geweke.diag(dep2.hh.mcmc);geweke.diag(dep4.hh.mcmc);geweke.diag(dep5.hh.mcmc);geweke.diag(dep7.hh.mcmc);geweke.diag(dep1.hh.lag.mcmc);geweke.diag(dep2.hh.lag.mcmc);geweke.diag(dep4.hh.lag.mcmc);geweke.diag(dep5.hh.lag.mcmc);geweke.diag(dep7.hh.lag.mcmc);geweke.diag(dep1.p.mcmc);geweke.diag(dep2.p.mcmc);geweke.diag(dep4.p.mcmc);geweke.diag(dep5.p.mcmc);geweke.diag(dep7.p.mcmc);geweke.diag(dep1.p.lag.mcmc);geweke.diag(dep2.p.lag.mcmc);geweke.diag(dep4.p.lag.mcmc);geweke.diag(dep5.p.lag.mcmc);geweke.diag(dep7.p.lag.mcmc)
#geweke.diag(dep1.hh.mcmc);geweke.diag(dep2.hh.mcmc);geweke.diag(dep3.hh.mcmc);geweke.diag(dep4.hh.mcmc);geweke.diag(dep5.hh.mcmc);geweke.diag(dep6.hh.mcmc);geweke.diag(dep7.hh.mcmc);geweke.diag(dep1.hh.lag.mcmc);geweke.diag(dep2.hh.lag.mcmc);geweke.diag(dep3.hh.lag.mcmc);geweke.diag(dep4.hh.lag.mcmc);geweke.diag(dep5.hh.lag.mcmc);geweke.diag(dep6.hh.lag.mcmc);geweke.diag(dep7.hh.lag.mcmc);geweke.diag(dep1.p.mcmc);geweke.diag(dep2.p.mcmc);geweke.diag(dep3.p.mcmc);geweke.diag(dep4.p.mcmc);geweke.diag(dep5.p.mcmc);geweke.diag(dep6.p.mcmc);geweke.diag(dep7.p.mcmc);geweke.diag(dep1.p.lag.mcmc);geweke.diag(dep2.p.lag.mcmc);geweke.diag(dep3.p.lag.mcmc);geweke.diag(dep4.p.lag.mcmc);geweke.diag(dep5.p.lag.mcmc);geweke.diag(dep6.p.lag.mcmc);geweke.diag(dep7.p.lag.mcmc)

###OUTPUT TABLES###
dep7.output<-cbind(
summary(dep7.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep7.hh.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep7.p.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep7.p.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)]
)
rownames(dep7.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep7.output,digits=4,caption="Models of dep7 from 1976-2008 with Four Different Gini Measures, Posterior Summaries from MCMC",label='dep7.table',align="lrrrrrrrrrrrr")

dep1.output<-cbind(
summary(dep1.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep1.hh.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep1.p.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep1.p.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)]
)
rownames(dep1.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep1.output,digits=4,caption="Models of dep1 from 1976-2008 with Four Different Gini Measures, Posterior Summaries from MCMC",label='dep1.table',align="lrrrrrrrrrrrr")

dep2.output<-cbind(
summary(dep2.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep2.hh.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep2.p.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep2.p.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)]
)
rownames(dep2.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep2.output,digits=4,caption="Models of dep2 from 1976-2008 with Four Different Gini Measures, Posterior Summaries from MCMC",label='dep2.table',align="lrrrrrrrrrrrr")

dep4.output<-cbind(
summary(dep4.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep4.hh.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep4.p.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep4.p.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)]
)
rownames(dep4.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep4.output,digits=4,caption="Models of dep4 from 1976-2008 with Four Different Gini Measures, Posterior Summaries from MCMC",label='dep4.table',align="lrrrrrrrrrrrr")

dep5.output<-cbind(
summary(dep5.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep5.hh.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep5.p.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)],
summary(dep5.p.lag.mcmc,quantiles=c(.05,.5,.95))$quantiles[2:11,c(2,1,3)]
)
rownames(dep5.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep5.output,digits=4,caption="Models of dep5 from 1976-2008 with Four Different Gini Measures, Posterior Summaries from MCMC",label='dep5.table',align="lrrrrrrrrrrrr")


#####ROBUSTNESS CHECKS FOR REVIEWERS#####

###BASE MODEL FOR NEW SPECIFICATIONS###
dep7.hh.mcmc<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.hh.mcmc,quantiles=c(.05,.5,.95))

###WITH PUBLIC OPINION###
dep7.hh.opinion<-MCMCregress(rAdjwelftotexptot~opinion+pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.hh.opinion,quantiles=c(.05,.5,.95))
#This variable is coded so that 1 means ``Government see to job and good standard of living" and 7 means ``Government let each person get ahead on his own". Hence, we would expect a negative correlation such that a higher average value (indicating a preference that government not get involved) would imply a more active government policy. However, we instead see a positive coefficient that does not have a robust effect.

opinion.output<-summary(dep7.hh.opinion,quantiles=c(.05,.5,.95))$quantiles[c(2:12),c(2,1,3)]
rownames(opinion.output)<-c("Public opinion","Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(opinion.output,digits=4,caption="Model of dep7 from 1976-2008 with Opinion",label='opinion.mod',align="lrrr")

###YEAR FIXED EFFECTS###
dep7.hh.year.fixed<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(year),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.hh.year.fixed,quantiles=c(.05,.5,.95))
year.output<-summary(dep7.hh.year.fixed,quantiles=c(.05,.5,.95))$quantiles[c(2:11),c(2,1,3)]
rownames(year.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(year.output,digits=4,caption="Model of dep7 from 1976-2008 with Year Fixed Effects",label='year.mod',align="lrrr")
year.output.side<-summary(dep7.hh.year.fixed,quantiles=c(.11,.5,.95))$quantiles[c(2),c(2,1,3)]; year.output.side #89%

###ALTERNATE DENOMINATOR###
#cor(fiscal$dv.gsp,fiscal$rAdjwelftotexptot,use="complete.obs")
#plot(y=fiscal$dv.gsp,x=fiscal$rAdjwelftotexptot)
dv.gsp.hh.mcmc<-MCMCregress(dv.gsp~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dv.gsp.hh.mcmc,quantiles=c(.05,.5,.95))
dv.gsp.output<-summary(dv.gsp.hh.mcmc,quantiles=c(.05,.5,.95))$quantiles[c(2:11),c(2,1,3)]
rownames(dv.gsp.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dv.gsp.output,digits=4,caption="Model Of the Ratio of Welfare Spending to GSP",label='gsp.denom',align="lrrr")

###ALTERNATIVE SPECIFICATION SUGGESTED BY REVIEWER THAT INCLUDES STATE DEBT###
debt.model<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+lnDebt+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(debt.model,quantiles=c(.05,.5,.95));geweke.diag(debt.model)

debt.output<-summary(debt.model,quantiles=c(.05,.5,.95))$quantiles[2:12,c(2,1,3)]
rownames(debt.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)","Log Total State Debt")
xtable(debt.output,digits=4,caption="Model of dep7 from 1976-2008 Including State Debt",label='table.debt',align="lrrr")

###CONSIDERING WHETHER THE EFFECT CHANGES OVER TIME###
dep7.hh.lm<-lm(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal);summary(dep7.hh.lm)#0.8350
alt.mod<-lm(rAdjwelftotexptot~pregini_hh_con+pregini_hh_con:as.numeric(year>=1986)+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal);summary(alt.mod)#0.8353
anova(dep7.hh.lm,alt.mod)

change.1986<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+pregini_hh_con:as.numeric(year>=1986)+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(change.1986,quantiles=c(.05,.5,.95));geweke.diag(change.1986)

change.output<-summary(change.1986,quantiles=c(.05,.5,.95))$quantiles[c(2:11,60),c(2,1,3)]
rownames(change.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)","Gini after 1986")
xtable(change.output,digits=4,caption="Model of dep7 from 1976-2008 with Time Switch",label='time.switch',align="lrrr")

###ALTERNATIVE SPECIFICATION THAT INCLUDES LAGGED PER CAPITA INCOME###
income.model<-MCMCregress(rAdjwelftotexptot~pregini_hh_con+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+lnpersincpc+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(income.model,quantiles=c(.05,.5,.95));geweke.diag(income.model)

income.output<-summary(income.model,quantiles=c(.05,.5,.95))$quantiles[2:12,c(2,1,3)]
rownames(income.output)<-c("Gini coefficient","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)","Log per capita income")
xtable(income.output,digits=4,caption="Model of dep7 from 1976-2008 Including Logged Per Capita Income",label='table.income',align="lrrr")

###ALTERNATIVE SPECIFICATION THAT USES 90% VS. MEDIAN###
dep7.topPerc<-MCMCregress(rAdjwelftotexptot~ratio9050+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.topPerc,quantiles=c(.05,.5,.95))
dep7.topPerc.output<-summary(dep7.topPerc,quantiles=c(.05,.5,.95))$quantiles[c(2:11),c(2,1,3)]
rownames(dep7.topPerc.output)<-c("Ratio of 90p to 50p","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep7.topPerc.output,digits=4,caption="Model of dep7 from 1976-2008 Using 90p to 50p Comparison",label='ratio9050',align="lrrr")

###ALTERNATIVE SPECIFICATION THAT USES MEAN VS. MEDIAN###
fiscal$meanMedian<-log(fiscal$AvgInc-fiscal$pre50p_hh_con)
dep7.meanMed<-MCMCregress(rAdjwelftotexptot~meanMedianRatio+gspgrowth+unemploy+lnpop+pop65+nonwhiteprop+gubelection+inst6010_nom+folded_ranney_4yrs+pre10p_hh_con+as.factor(statefip),data=fiscal,burnin=15000,mcmc=150000);summary(dep7.meanMed,quantiles=c(.05,.5,.95))
dep7.meanMed.output<-summary(dep7.meanMed,quantiles=c(.05,.5,.95))$quantiles[c(2:11),c(2,1,3)]
rownames(dep7.meanMed.output)<-c("Ratio of mean to median","Change in GSP","Unemployment","Log population","Pop. over 65","Nonwhite pop.","Gubernatorial elec.","Government ideology","Party competition","10p income (1000s)")
xtable(dep7.meanMed.output,digits=4,caption="Model of dep7 from 1976-2008 Using Mean vs.~Median Comparison",label='meanMedian',align="lrrr")

###POSTERIOR DENSITY PLOTS###
#pdf('density7.pdf',width=3,height=3,pointsize=7)
plot(density(dep7.hh.mcmc[,2]),xlab="Effect of Gini Coefficient",ylab="Density",main="",xlim=c(-.1,.7),ylim=c(0,14),lty=1)
lines(density(dep7.hh.lag.mcmc[,2]),col='red',lty=2)
lines(density(dep7.p.mcmc[,2]),col='blue',lty=3)
lines(density(dep7.p.lag.mcmc[,2]),col='darkgreen',lty=4)
legend(x=-.1,y=14,col=c("black","red","blue","darkgreen"),lty=c(1,2,3,4),legend=c("Household","Lag household","Personal","Lag personal"),ncol=2)
#dev.off()

#pdf('density1.pdf',width=3,height=3,pointsize=7)
plot(density(dep1.hh.mcmc[,2]),xlab="Effect of Gini Coefficient",ylab="Density",main="",xlim=c(-1,9.5),ylim=c(0,2.25),lty=1)
lines(density(dep1.hh.lag.mcmc[,2]),col='red',lty=2)
lines(density(dep1.p.mcmc[,2]),col='blue',lty=3)
lines(density(dep1.p.lag.mcmc[,2]),col='darkgreen',lty=4)
legend(x=-.05,y=2.25,col=c("black","red","blue","darkgreen"),lty=c(1,2,3,4),legend=c("Household","Lag household","Personal","Lag personal"),ncol=2)
#dev.off()

#pdf('density2.pdf',width=3,height=3,pointsize=7)
plot(density(dep2.hh.mcmc[,2]),xlab="Effect of Gini Coefficient",ylab="Density",main="",xlim=c(-1,9.5),ylim=c(0,2.25),lty=1)
lines(density(dep2.hh.lag.mcmc[,2]),col='red',lty=2)
lines(density(dep2.p.mcmc[,2]),col='blue',lty=3)
lines(density(dep2.p.lag.mcmc[,2]),col='darkgreen',lty=4)
legend(x=-.05,y=2.2,col=c("black","red","blue","darkgreen"),lty=c(1,2,3,4),legend=c("Household","Lag household","Personal","Lag personal"),ncol=2)
#dev.off()

#pdf('density4.pdf',width=3,height=3,pointsize=7)
plot(density(dep4.hh.mcmc[,2]),xlab="Effect of Gini Coefficient",ylab="Density",main="",xlim=c(-.1,.6),ylim=c(0,16),lty=1)
lines(density(dep4.hh.lag.mcmc[,2]),col='red',lty=2)
lines(density(dep4.p.mcmc[,2]),col='blue',lty=3)
lines(density(dep4.p.lag.mcmc[,2]),col='darkgreen',lty=4)
legend(x=-.1,y=16.2,col=c("black","red","blue","darkgreen"),lty=c(1,2,3,4),legend=c("Household","Lag household","Personal","Lag personal"),ncol=2)
#dev.off()

#pdf('density5.pdf',width=3,height=3,pointsize=7)
plot(density(dep5.hh.mcmc[,2]),xlab="Effect of Gini Coefficient",ylab="Density",main="",xlim=c(-.1,.7),ylim=c(0,14),lty=1)
lines(density(dep5.hh.lag.mcmc[,2]),col='red',lty=2)
lines(density(dep5.p.mcmc[,2]),col='blue',lty=3)
lines(density(dep5.p.lag.mcmc[,2]),col='darkgreen',lty=4)
legend(x=-.1,y=14,col=c("black","red","blue","darkgreen"),lty=c(1,2,3,4),legend=c("Household","Lag household","Personal","Lag personal"),ncol=2)
#dev.off()

###INTERPRETATION###
#standard deviation rise
0.5019*sd(fiscal$pregini_hh_con,na.rm=T)

#10th to 90th percentile effect
0.5019*(0.5308072-0.4338348)

#Using dep1, how much wealth redistribution would we see?
summary(fiscal$lnWelfTotPC[fiscal$year==2007])
summary(fiscal$lnpersincpc[fiscal$year==2007])
summary(fiscal$pre10p_hh_con[fiscal$year==2007])
summary(fiscal$pre50p_hh_con[fiscal$year==2007])

#increase<-7.7975*(0.5992-0.3789); increase #percentage increase in spending for min to max Gini
increase<-7.7975*(0.5308072-0.4338348); increase #percentage increase in spending for 10p to 90p Gini
exp(7.099)*(increase) #dollar increase at the average level of spending
exp(7.099) #average level of spending
exp(7.099)*(1+increase) #level of spending if Gini rose

exp(10.55) #average per capita income
exp(7.099)*(increase)/exp(10.55) #as a percentage of average per capita income
exp(7.099)*(increase)/793.90 #as a percentage of average 10th percentile income
